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•  Stochastic,  3D  carbon  fiber  material  models  generated  with  controlled  porosity  heterogeneity. 

•  Lattice  Boltzmann  method  employed  to  measure  bulk  hydrodynamic  properties. 

•  Porosity  heterogeneity  contributed  a  negligible  effect  on  modeled  transport,  compared  to  porosity  and  material  composition. 

•  Tortuosity  increases  with  binder  &  PTFE  fraction,  which  was  unexpected  by  the  authors. 

•  Permeability  increases  with  permeability  and  binder  &  PTFE  fraction. 
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In  this  paper,  we  study  the  effect  of  porosity  heterogeneity  on  the  bulk  hydrodynamic  properties 
(permeability  and  tortuosity)  of  simulated  gas  diffusion  layers  (GDLs).  The  porosity  distributions  of  the 
heterogeneous  reconstructed  samples  are  similar  to  those  previously  reported  in  the  literature  for  Toray 
TGP-H  120™  GDLs.  We  use  the  lattice  Boltzmann  method  to  perform  pore-level  flow  simulations  in  the 
reconstructed  GDL  samples.  Using  the  results  of  pore-level  simulations,  the  effect  of  porosity  distribution 
is  characterized  on  the  predicted  in-  and  cross-plane  permeability  and  tortuosity.  It  was  found  that 
porosity  heterogeneity  causes  a  higher  in-plane  permeability  and  lower  in-plane  tortuosity,  while  the 
effect  is  opposite  in  the  cross-plane  direction,  that  is  a  lower  cross-plane  permeability  and  a  higher  cross¬ 
plane  tortuosity. 

We  further  investigate  the  effect  of  adding  poly-tetra-fluoro-ethylene  (PTFE)  &  binder  material  to  the 
reconstructed  GDL  samples.  Three  fiber  volume  percentages  of  50,  75,  and  100%  are  considered.  Overall, 
increasing  the  fiber  volume  percentage  reduces  the  predicted  in-  and  cross-plane  permeability  and 
tortuosity  values.  A  previously  reported  relationship  for  permeability  of  fibrous  materials  is  fitted  to  the 
predicted  permeability  values,  and  the  magnitude  of  the  fitting  parameter  is  reported  as  a  function  of 
fiber  volume  percentage. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Hydrogen-based  power  systems,  specifically  fuel  cells,  have 
attracted  significant  attention  in  the  last  decade  due  to  their  low 
local  CO2  emission,  high  efficiency,  and  versatility.  Among  different 
types  of  fuel  cells,  polymer  electrolyte  membrane  (PEM)  fuel  cells, 
which  typically  have  energy  conversion  efficiencies  over  50%,  are  a 
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promising  option  for  the  transportation  sector,  distributed  power 
generation,  and  backup  power  generation,  to  name  a  few  applica¬ 
tions  [1-4], 

Considerable  advancements  have  been  made  towards  the 
commercial  development  of  PEM  fuel  cells  in  areas  such  as  reduced 
platinum  loading,  higher  current  density,  and  reduced  costs; 
however,  there  still  exist  challenges  that  should  be  removed  before 
PEM  fuel  cells  can  become  commercially  viable.  Mass  transport 
issues  such  as  liquid  water  flooding  at  high  current  working  con¬ 
ditions  leads  to  inefficient  performance  [5],  Overcoming  these 
operational  challenges  requires  more  research  on  characterization 
and  performance-evaluation  of  PEM  fuel  cell  materials,  which  will 
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inform  new  material  designs  for  improved  efficiency  and  product 
reliability. 

Computational  modeling,  next  to  experimental  studies,  has 
been  widely  used  to  study  mass,  energy,  and  charge  transport  at 
different  length  scales  inside  PEM  fuel  cells.  The  numerical  studies 
can  be  categorized  into  two  groups,  namely:  (i)  studies  that  use  first 
principles,  thus,  can  only  include  small  system  sizes  comparable  to 
length  scales  of  the  transport  phenomenon  under  study,  and  (ii) 
studies  that  use  models  and  relations  to  predict  transport  phe¬ 
nomena  inside  multiple  components  of  a  fuel  cell,  and  predict  the 
cell-level  behavior.  These  models  usually  include  larger  combined 
domains,  and  need  effective  transport  properties  of  each  part  of  the 
cell  as  an  input.  The  aim  of  many  of  the  works  under  the  first 
category  is  to  predict  effective  transport  properties  of  individual 
components  of  a  PEM  fuel  cell,  which  will  be  used  in  the  studies  of 
the  second  type. 

One  particular  part  of  a  PEM  fuel  cell,  which  is  the  focus  of  this 
paper,  is  the  gas  diffusion  layer  (GDL).  The  GDL  is  a  layer  consisting 
of  carbon  fibers  in  either  woven  or  non-woven  formats.  The 
thickness  of  a  GDL  is  typically  between  200—400  pm,  with  fiber 
diameters  in  the  range  of  7— 10  pm.  In  a  conventional  PEM  fuel  cell, 
the  GDL  is  placed  between  the  micro  porous  layer  (MPL)  and  a 
graphite  plate  that  includes  gas  channels.  The  produced  water  in¬ 
side  the  cell  passes  through  the  GDL  to  be  guided  out  of  the  system 
via  gas  channels.  The  GDL  should  be  of  overall  high  thermal  and 
electrical  conductivity,  and  should  be  structured  in  a  way  that 
prevents  excess  water  accumulation. 

The  GDL  influences  many  aspects  of  a  PEM  fuel  cell’s  operation, 
thus,  accurate  prediction  of  its  effective  transport  properties  is 
crucial  in  understanding  the  cell  performance.  El-Kharouf  et  al.  [6] 
presented  a  comprehensive  review  of  available  ex-situ  measure¬ 
ment  techniques  for  determining  a  range  of  GDL  properties. 

A  review  of  available  methods  for  bulk  porosity  measurement 
techniques  is  presented  by  Cindrella  et  al.  [7],  Determining  effective 
thermal  conductivity  of  GDL  samples  has  also  been  subject  of  many 
previous  studies  [8  3).  Another  important  bulk  transport  prop¬ 

erty  of  a  GDL  is  its  absolute  permeability.  Bulk  permeability  pre¬ 
diction  in  GDL  samples  has  been  widely  investigated  numerically 
and  experimentally  [14—19],  All  these  works  focus  on  determining 
the  bulk  permeability,  and  do  not  discuss  the  effect  of  pore-level 
structural  parameters  on  the  permeability.  Recently,  Fishman 
et  al.  [20]  showed  that  porosity  distribution  through  the  thickness 
of  GDL  samples  is  not  uniform.  Through  rigorous  analysis  of 
computed  tomography  (CT)  scan  images  of  samples  from  different 
manufacturers,  they  reported  typical  porosity  distribution  profiles 
of  different  types  of  GDLs.  In  another  work,  Fishman  and  Bazylak 
[21]  studied  the  effect  of  poly-tetra-fluoro-ethylene  (PTFE)  and 
binder  material  addition  on  the  overall  through-thickness  porosity 
profiles  of  GDL  samples.  They  reported  that  PTFE  and  binder  ma¬ 
terial  accumulate  in  low  porosity  areas,  where  the  density  of  fibers 
is  higher. 

It  is  known  that  the  heterogeneous  porosity  distribution  affects 
the  overall  transport  properties  of  porous  material  [22,23],  how¬ 
ever,  this  dependance  has  not  been  quantified.  Fishman  et  al.  [24] 
applied  the  available  relations  for  permeability/tortuosity  of 
porous  materials  to  the  measured  porosity  profile  inside  the  GDL, 
and  mapped  the  porosity  profile  into  permeability/tortuosity  pro¬ 
files.  In  this  work,  we  show  that  their  method  of  calculating  the 
resultant  bulk  permeability  is  only  accurate  for  the  in-plane 
direction. 

In  this  work,  we  use  pore-level  flow  simulations  to  investigate 
how  the  heterogeneous  porosity  distribution  of  the  GDL  affect 
overall  bulk  hydrodynamic  properties  (i.e.,  permeability  and  tor¬ 
tuosity).  In  order  to  conduct  a  parametric  study,  GDL  samples  are 
reconstructed  over  the  entire  range  of  relevant  bulk  porosity  values 


with  homogeneous  and  heterogeneous  porosity  profiles.  The  het¬ 
erogeneous  profile  is  similar  to  that  reported  previously  by  Fish¬ 
man  et  al.  [20,21]  from  CT  scan  imaging  of  core  section  of  Toray 
TGP-H  120™  GDL  samples.  Furthermore,  the  effect  of  the  PTFE  & 
binder  addition  on  bulk  hydrodynamic  properties  of  GDL  samples 
has  also  been  systematically  investigated.  PTFE  and  binder  are 
treated  as  one  solid  material. 

We  have  followed  the  lattice  Boltzmann  method  (LBM)  [25] 
approach  for  solving  the  steady-state  three-dimensional  single¬ 
phase  fluid  flow  in  the  reconstructed  geometries.  Due  to  its 
capability  in  easily  handling  complex  geometries,  the  LBM  has 
been  used  widely  to  model  fluid  flow  in  reconstructed  porous 
structures  [17,26-29],  In  a  recent  work,  Froning  et  al.  used  the 
same  methodology  to  study  the  compression  effect  of  the  pre¬ 
dicted  permeability  in  three-dimensional  reconstructed  media 
[30],  In  all  these  works,  the  porous  structure  is  reconstructed  by 
randomly  placing  cylindrical  fibers  in  a  three-dimensional 
domain.  However,  while  randomly  placing  the  fibers  in  the 
planes,  we  control  the  through-plane  porosity  profile.  This  enables 
us  to  systematically  study  the  effect  of  porosity  heterogeneity  on 
bulk  properties.  Furthermore,  we  also  investigate  the  effect  of 
Binder  &  PTFE  addition. 

2.  Methodology 

In  this  work,  steady-state  three-dimensional  fluid  flow  was 
simulated  using  LBFlow,1  which  is  a  single-relaxation-time  imple¬ 
mentation  of  the  LBM  [25]. 

The  LBM  is  a  mesoscopic  approach  for  modeling  fluid  hydro¬ 
dynamics,  where  fictitious  particles  represent  fluid  volumes.  The 
macroscopic  hydrodynamic  parameters  (e.g.,  density  and  velocity) 
are  calculated  based  on  the  moments  of  the  distribution  function 
of  these  fictitious  particles  on  each  lattice  site.  The  LBM  algorithm 
includes  two  steps,  namely:  (i)  streaming  step,  where  the  distri¬ 
bution  functions  move  in  the  direction  of  their  assigned  discrete 
velocity  set,  and  (ii)  the  collision  step,  where  new  distributions  are 
calculated  on  each  lattice  site  based  on  pre-defined  collision  rules. 
The  collision  step  includes  the  effect  of  the  external  forces  and  the 
boundary  conditions.  The  details  of  the  numerical  methodology 
and  method  validation  for  LBFlow  are  presented  elsewhere 
[27,31-33], 

The  computational  domain  of  the  porous  region  is  specified  as  a 
three-dimensional  binary  array  of  solid  and  fluid  nodes.  Physical 
properties  of  working  fluid  are  set  to  those  of  water  at  20  °C.  The 
coordinate  system  is  defined  such  that  the  z-direction  is  aligned 
with  the  cross-plane  direction,  and  the  x—y  plane  is  the  plane 
parallel  to  the  carbon  fibers  of  the  GDLs.  A  pressure  gradient,  in  the 
form  of  a  body  force,  is  applied  in  the  direction  of  each  main  axis, 
and  the  obtained  three-dimensional  flow  field  is  averaged  to  yield  a 
volume-averaged  velocity  vector,  ii.  The  permeability  tensor  (K)  is 
then  calculated  using  Darcy’s  law,  as  follows: 


where  u,  is  the  volume-averaged  velocity  in  the  direction  i,  Ky  is  the 
ijth  component  of  the  permeability  tensor,  ji  is  the  working  fluid 
viscosity,  and  dp/dxj  is  the  applied  pressure  gradient  in  the  direc¬ 
tion  j.  Considering  the  structure  of  GDL  samples  in  this  work,  we 
assumed  that  off  diagonal  elements  of  the  permeability  tensor  (% 
I  ¥=  j )  has  negligible  value  in  comparison  to  its  diagonal  elements 


1  Available  from  http://www.lbflow.co.uk. 
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Another  important  bulk  property  of  porous  materials  is  the 
tortuosity  (t).  Tortuosity  is  usually  related  to  the  effective  diffusion 
coefficient  (Detr)  of  porous  media  as  Deg  =  0/tD,  where  D  is  the 
binary  diffusion  coefficient  [7],  We  used  the  following  formula  to 
calculate  the  tortuosity  in  an  arbitrary  direction  i  [27]: 


Ti  E£™S 


(2) 


where  umag  is  the  velocity  magnitude  calculated  as  +  uj  +  u^, 
m  is  the  velocity  component  in  direction  i,  and  the  summation  is 
performed  over  all  the  fluid  nodes. 


3.  Geometry  reconstruction 

To  simulate  fluid  and  mass  transport  in  a  GDL-like  material, 
digital  models  are  created  using  stochastic  modeling  techniques  for 
fibrous  materials.  An  Nx  x  Ny  x  Nz  voxelized  domain  is  filled  with 
cylinders  until  a  set  desired  porosity  is  reached.  The  sample  gen¬ 
eration  technique  that  we  used  in  this  work  most  closely  resembles 
that  previously  reported  in  Refs.  [34,35],  but  follows  similar 
modeling  steps  to  the  family  of  methods  previously  applied  for 
digital  GDL  reconstructions  [12,27,36-39],  Similar  to  [12,39], 
samples  are  assumed  to  be  primarily  composed  of  straight,  cylin¬ 
drical  fibers,  oriented  along  the  primary  (x—y)  plane  of  the  material, 
where  fibers  are  allowed  to  intersect.  Additionally,  to  create  het¬ 
erogeneous  samples,  the  cross-plane  porosity  distribution  of  the 
core  section  of  Toray  TGP-H  120™  with  10%  wt  PTFE  [21  ]  is  applied 
to  a  group  of  samples,  following  the  methods  described  in  Ref.  [40], 
The  cross-plane,  local  porosity  distribution  of  a  single  recon¬ 
structed  material  is  displayed  in  Fig.  1,  and  compared  with  the 
distribution  used  as  input  for  material  generation.  The  sharp  fea¬ 
tures  in  generated  material’s  profile  are  a  result  of  the  0.42  mm2 
sample  size,  compared  to  the  25  mm2  material  scanned. 

We  also  created  samples  with  added  binder  and  PTFE.  The 
combination  of  these  two  coatings  is  represented  as  a  secondary 
material,  residing  in  the  tight  crevices  near  fiber  intersections, 
following  the  methods  of  [36,39],  The  fiber  volume  percentage,  a,  is 
calculated  as  the  ratio  of  the  fibers’  volume  to  the  total  volume  of 
the  solid  phase. 

Fig.  2  displays  three  reconstructed  samples.  Each  sample  is 
made  to  match  the  same  specified  cross-plane  porosity  distribution 
and  have  an  average  bulk  porosity  of  0.70.  However,  the  fiber  to 
binder  &  PTFE  ratio  is  varied  between  the  three  samples.  For  a 
better  comparison,  the  three  materials  were  generated  from  the 
same  list  of  fiber  positions  and  orientations,  however,  those  with 


Cross-Plane  Position  (urn) 


Fig.  1.  Comparison  between  the  cross-plane  porosity  distribution  used  as  input 
(dashed)  to  the  material  generation,  and  that  of  a  single,  100%  fibrous  material 
generated  for  this  study  (solid). 


more  binder  &  PTFE  employed  a  truncated  list.  In  Fig.  2(a),  the  fiber 
volume  percentage  is  50%.  In  Fig.  2(b),  the  solid  matrix  is  75%  fibers. 
Finally,  in  Fig.  2(c),  there  is  no  binder  &  PTFE  at  all.  Fig.  2(d-f) 
provide  a  cross-sectional  image  of  the  reconstructed  samples.  At 
this  cross-plane  position,  the  local  porosity  is  around  0.65.  As  can 
be  seen  from  these  cross-sections,  the  inclusion  of  more  fibers 
creates  a  higher  number  of  smaller  pores.  The  accumulation  of 
binder  &  PTFE  (marked  in  blue  (in  the  web  version))  at  the  tight 
crevices  near  fiber  intersections  is  visible  in  Fig.  8(d)  and  (e). 

Two  considerations  that  were  made  specifically  for  the  sample 
reconstructions  in  this  work  should  be  highlighted.  First,  the  fiber 
placement  was  done  in  such  a  way  as  to  provide  ideal  sidewall 
boundaries  for  periodic  simulations.  This  means  that  when  a 
sidewall  is  mated  with  its  opposite  face,  the  boundaries  are  indis¬ 
tinguishable  from  the  internal  volume  of  the  material.  To  accom¬ 
plish  this  feature,  individual  fiber  placement  was  accomplished  in  4 
steps: 

1 )  A  cylindrical  fiber  of  length  Nx  is  placed  in  the  center  of  the 
domain,  such  that  its  center  point  is  at  JVx/2,  JVy/2,  Nz/ 2, 
extending  only  in  the  x-direction. 

2)  The  fiber  is  translated  vertically  a  random  amount,  weighted  by 
the  intended  porosity  distribution,  as  described  in  Ref.  [34], 

3)  The  fiber  is  then  rotated  in  the  x—y  plane  at  a  random  angle 
between  0  and  180,  about  its  center  point. 

4)  The  fiber  is  then  translated  again  by  a  random  distance  in  both 
the  x  and  y  directions,  such  that  the  center  is  moved  to  another 
point  in  the  x—y  plane,  with  a  uniform  probability. 

5)  Any  portion  of  the  fiber  that  extends  past  the  edge  of  the 
domain,  is  retained,  and  now  extends  into  the  opposite  face  of 
the  domain,  such  that  if  you  matched  opposite  faces  together,  a 
straight  fiber  of  length  of  Nx  would  be  formed. 

The  size  of  the  reconstructed  domain  is  Nx  =  325  and  Ny  —  325 
lattice  units  in  x  andy  directions.  The  sample  thickness  has  changed 
based  on  different  models,  while  Nz  stays  around  150  lattice  units 
typically.  The  authors  have  previously  studied  the  effect  of  domain 
size  on  the  predicted  permeability  of  porous  structures  [23],  It  has 
been  shown  that  this  domain  size  is  large  enough  to  yield  lattice 
independent  results.  The  fibers  are  4  lattice  units  in  diameter.  The 
effect  of  fiber  diameter  on  the  predicted  permeability  has  also  been 
elaborated  in  a  previous  work  of  the  authors,  where  the  fiber  di¬ 
ameters  were  systematically  changed  and  predicted  permeability 
was  reported  [27], 

The  second  consideration  that  was  made  for  this  work  pertained 
to  the  addition  of  the  binder  &  PTFE  material  onto  the  placed  fibers 
such  that  a  desired  porosity  distribution  would  be  matched.  The 
porosity  distributions  used  in  this  work  are  cross-plane,  meaning 
that  the  local  porosity  of  the  generated  samples  should  be  a  func¬ 
tion  of  the  z  coordinate  alone.  Because  the  binder  &  PTFE  is  placed 
as  a  wetting  fluid,  it  disproportionally  resides  in  regions  of  low 
porosity.  Therefore,  although  this  scheme  provides  global  control  of 
how  much  binder  &  PTFE  gets  added  to  the  sample,  there  is  no  local 
control  for  binder  &  PTFE  placement.  To  remedy  this,  the  z-position 
weighting  function  (J[z))  for  random  fiber  placement  is  modified. 
The  original  function  of: 

f(z)  =  1  -  *(z)  (3) 

where  </>  is  the  porosity,  assumes  all  material  in  the  sample  will  be 
placed  in  that  fashion.  We  investigated  a  general  function,  /,  that 
could  provide  an  appropriate  weighting  for  cross-plane  fiber 
placement,  while  accounting  for  the  fiber  to  binder  &  PTFE  ratio,  a, 
as  well.  As  a  result  of  the  investigation,  it  was  apparent  that  not  only 
a,  but  also  the  shape  of  the  original  function,/,  played  an  important 
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(a)  50%  fibre  (b)  75%  fibre  (c)  100%  fibre 


(d)  50%  fibre  (e)  75%  fibre  (f)  100%  fibre 


Fig.  2.  (Top)  three-dimensional  structure  and  (bottom)  two-dimensional  slices  of  reconstructed  GDL  samples  with  three  different  fiber  percentages  of  50%,  75%,  and  100%.  The  color 
blue  shows  the  added  binder  &  PTFE  material,  while  the  white  color  is  representative  of  fibers.  Average  bulk  porosity  of  all  three  samples  is  0.70,  while  the  local  porosity  of  the 
shown  slices  is  around  0.65.  (a)  50%  fiber;  (b)  75%  fiber;  (c)  100%  fiber;  (d)  50%  fiber;  (e)  75%  fiber;  (f)  100%  fiber.  (For  interpretation  of  the  references  to  color  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  article.) 


role.  Although  no  general  function  was  found,  through  trial  and 
error,  a  relatively  good  approximation  of  /  was  arrived  at  using  the 
following  function: 

m=nz)(^\f^  +i)  (4) 


4.  Results 

4.3.  Introduction 

We  used  a  pore-level  flow  modeling  approach  to  study  the  effect 
of  simulated  GDLs’  structure  on  their  hydrodynamic  bulk  proper¬ 
ties,  namely:  permeability  and  tortuosity.  Two  main  structural 
features  of  GDL  samples  are  investigated  in  detail.  First,  effect  of 
porosity  heterogeneity  on  the  overall  permeability  and  tortuosity  of 
samples  is  studied.  Second,  the  effect  of  the  addition  of  binder  & 
PTFE  material  on  the  bulk  properties  is  studied  for  different  values 
of  fiber  percentage. 

4.2.  Effect  of  porosity  heterogeneity  in  fiber-only  samples 

In  order  to  study  the  effect  of  porosity  distribution  on  the  pre¬ 
dicted  permeability  and  tortuosity  of  reconstructed  GDL  samples, 
we  created  samples  of  homogeneous  and  heterogeneous  porosity 
distributions  in  the  through-plane  direction  for  the  entire  range  of 
relevant  porosity,  i.e.,  bulk  porosities  between  0.50  and  0.95.  For  the 
heterogeneous  samples,  we  chose  the  core  region  of  the  porosity 
profile  of  Toray  TGP-H  120™  with  10%  wt.  PTFE  samples,  as 


reported  in  Ref.  [21  ].  This  profile  was  chosen  as  it  provides  the  most 
accentuated  porosity  heterogeneity,  among  all  the  samples  that 
were  reported  in  Refs.  [20,21  ].  At  this  stage,  all  the  samples  consist 
of  only  fibers  and  no  secondary  material  (i.e.,  binder  &  PTFE)  is 
added. 

The  normalized  predicted  values  of  in-  and  cross-plane 
permeability  for  homogeneous  and  heterogeneous  GDL  samples 
are  presented  in  Fig.  3.  The  permeability  values  are  normalized  by 
the  square  of  fibers’  radius.  It  can  be  seen  that  for  the  in-plane 
permeability,  porosity  heterogeneity  increases  the  overall  perme¬ 
ability  of  the  sample:  consistent  with  what  previously  discussed  in 
Ref.  [22],  However,  for  the  cross-plane  direction,  the  porosity  het¬ 
erogeneity  decreases  the  overall  permeability  of  the  samples.  One 


Porosity,  <)> 

Fig.  3.  Predicted  values  of  normalized  in-  and  cross-plane  permeability  of  homoge¬ 
neous  and  heterogeneous  reconstructed  GDL-like  samples  as  a  function  of  porosity. 
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and  cross-plane  permeability  of  heterogeneous  GDL-like  samples. 


should  not  that,  although  the  effect  of  porosity  heterogeneity  is 
consistent  for  the  entire  range  of  porosity,  it  is  however  negligible 
compared  to  the  effect  of  other  parameters,  such  as  porosity  and 
fiber  aspect  ratio  [27], 

In  an  attempt  to  analytically  predict  the  effective  permeability  of 
heterogeneous  samples  from  their  porosity  profile,  we  assume  the 
GDL  is  made  up  of  parallel  layers,  as  schematically  shown  in  Fig.  4. 
These  imaginary  planes  are  supposed  to  be  homogeneous,  with  a 
uniform  porosity  corresponding  to  their  respective  location  on  the 
cross-plane  porosity  profile.  Knowing  the  porosity  of  each  layer,  its 
permeability  can  be  calculated  using  available  relations  in  the 
literature. 

For  each  layer,  one  can  write  the  Darcy’s  law  as  Qj  =  UjAj  [KjAP,/ 
/fLj]Aj  for  either  in-plane  or  cross-plane  directions.  Q,  and  u,  are  the 
flow  rate  and  average  velocity  through  layer  i  in  the  direction  of 
interest.  K,  and  AP,  are  the  permeability  and  pressure  drop,  both  in 
the  direction  of  the  mean  flow.  L*  is  the  length  of  the  layer  i  in  the 
direction  of  mean  flow,  and  ji  is  the  dynamic  viscosity  of  the 
working  fluid.  A,-  is  the  cross  section  of  the  layer  i  normal  to  the 
direction  of  interest.  It  should  be  noted  that  as  each  layer  is  sup¬ 
posed  to  be  homogeneous,  the  K,-  is  the  same  in  all  directions. 

For  the  in-plane  direction,  one  can  easily  observe  that 
APmt  =  AP,,  Qtot  =  £Qi,  Atot  =  J2a i,  and  hot  =  4  Therefore,  the 
overall  in-plane  permeability  can  be  obtained  as  a  function  of  the 
permeability  of  n  imaginary  layers  that  are  stacked  on  top  of  each 
other,  as  follows: 


=  (5) 

where  n  is  the  number  of  layers. 

For  the  cross-plane  direction,  flow  direction  is  normal  to  the 
planes,  and  thus,  one  can  see  that  APtot  =  ]TAPi,  Qtot  =  Q,,  Atot  =  A,-, 
and  Ltot  =  24  Therefore,  the  overall  cross-plane  permeability  is 
calculated  as  follows: 


l<cp 


n 

EUW 


(6) 


Table  1  presents  a  comparison  of  the  numerically  predicted 
values  of  in-  and  cross-plane  permeability  against  those  predicted 
from  Eqs.  (5)  and  (6)  for  two  representative  porosities  of  0.75  and 
0.85.  We  used  the  relationship  of  Nabovati  et  al.  [27]  to  predict  the 
permeability  of  each  layer  using  its  corresponding  porosity.  It  can 
be  seen  that  above-mentioned  relations  can  reasonable  predict  the 
permeability  of  heterogeneous  samples,  if  the  permeability  of  each 


Table  1 

Comparison  of  the  numerically  predicted  in-  and  cross-plane  permeability  values 
against  those  predicted  from  Eqs.  (5)  and  (6)  and  relationship  of  Nabovati  et  al.  [26] 
for  two  representative  porosities  of  0.71  and  0.80. 


Overall  Kip  Kcp 

porosity  Numerical  Analytical,  Numerical  Analytical, 

_ Eq.  (5) _ Eg,  (6) 

0.71  6.38  x  10“12  5.52  x  10  12  3.33  x  10  12  3.77  x  10  12 

0.80  16.11  x  10“12  15.04  x  10-12  8.96  x  10“12  9.88  x  1CT12 


imaginary  layer  is  correlated  to  its  porosity  using  available  relations 
in  the  literature. 

The  predicted  values  of  tortuosity  for  samples  of  homogeneous 
and  heterogeneous  porosity  distributions  are  shown  in  Fig.  5  as  a 
function  of  porosity.  It  can  be  seen  that  the  cross-plane  tortuosity 
values  are  consistently  larger  than  the  in-plane  tortuosity  values  for 
both  homogeneous  and  heterogeneous  porosity  distributions.  This 
behavior  is  opposite  to  that  of  permeability  values.  Furthermore,  it 
can  be  seen  that  the  heterogeneity  of  porosity  distribution  in¬ 
creases  the  cross-plane  tortuosity  and  decreases  the  in-plane  tor¬ 
tuosity.  However,  similar  to  the  predicted  values  of  permeability  in 
Fig.  3,  the  effect  of  porosity  distribution  on  the  predicted  values  of 
tortuosity  is  negligible  in  comparison  to  the  effect  of  fibers’  direc¬ 
tion  (cross-plane  vs.  in-plane)  and  porosity. 

4.3.  Effect  of  binder  &  PTFE 

In  this  section,  we  study  the  effect  of  the  addition  of  binder  & 
PTFE  to  the  base  fibrous  structure  of  GDLs  on  the  predicted 
permeability  and  tortuosity  values. 

Fig.  6  shows  the  normalized  predicted  values  of  in-  and  cross¬ 
plane  permeability  as  a  function  of  porosity  for  three  different  fi¬ 
ber  volume  percentages  of  50,  75,  and  100%.  It  can  be  seen  that  for 
both  directions,  decreasing  the  fiber  percentage  (adding  more 
binder  &  PTFE)  increases  the  permeability.  The  effect  of  fiber  per¬ 
centage  on  permeability  is  consistent  over  the  entire  range  of 
examined  porosity. 

The  positive  effect  of  adding  binder  &  PTFE  on  the  permeability 
can  be  explained  considering  the  difference  in  geometrical  shapes 
of  fibers  and  binder  &  PTFE  materials,  as  shown  in  Fig.  2.  Fibers  are 
in  the  form  of  straight  long  circular  cylinders;  while  the  added 
binder  &  PTFE  materials  are  more  in  the  form  of  bulk  of  materials 
accumulated  at  the  intersection  of  fibers  in  a  usually  irregular  bulk¬ 
like  form.  The  fibers,  thus,  have  a  higher  specific  surface  area 
(surface  area  per  unit  volume)  than  the  added  binder  &  PTFE 


Fig.  5.  Predicted  values  of  overall  in-  and  cross-plane  tortuosity  for  homogeneous  and 
heterogeneous  reconstructed  GDL-like  samples. 
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Fig.  6.  Predicted  values  of  normalized  in-plane  (top)  and  cross-plane  (bottom) 
permeability  of  reconstructed  GDL  samples  with  heterogeneous  porosity  distribution 
and  three  different  fiber  percentages  of  50%,  75%,  and  100%. 


material.  Samples  of  the  same  porosity  but  of  different  volume 
percentages  of  fibers  will  all  have  the  same  volume  of  materials; 
however,  because  of  different  specific  surface  area,  those  with  a 
higher  fiber  volume  percentage  (lower  binder  &  PTFE  volume 
percentage)  will  have  a  higher  solid  surface  area  that  is  exposed  to 
the  working  fluid.  The  reduction  in  the  permeability  is  mainly 
because  of  the  friction  between  the  solid  material  and  the  flowing 
fluid.  Thus,  a  higher  specific  surface  area,  means  a  larger  surface  for 
friction  to  act,  which  results  in  a  lower  permeability. 

To  quantify  the  effect  of  fiber  volume  percentage  on  the  pre¬ 
dicted  permeability  values,  we  modify  a  previously  reported  rela¬ 
tion  for  the  permeability  of  three-dimensional  random  fibrous 
structures  by  Nabovati  et  al.  [27]: 


Predicted  values  of  Cj  in  Eq.  (7)  for  in-  and  cross-plane  directions  for  different  values 
of  fiber  volume  percentage  from  curve  fitting  to  the  numerically  predicted  values  of 
permeability  for  samples  with  porosities  smaller  than  0.85.  0C  and  C2  are  constant 
and  equal  to  those  reported  by  Nabovati  et  al.  [26],  ie„  4 Sc  =  0.074  and  C2  =  2.310. 

In-plane  Cross-plane 

50%  fiber  volume  1.1480  0.6220 

75%  fiber  volume  0.8183  0.4625 

100%  fiber  volume  0.7057  0.3806 

Nabovati  et  al.  [26]  0.491 


where  a  is  the  fiber  radius,  (f>  is  porosity,  and  0C  is  the  percolation 
threshold.  Ci  and  C2  are  two  fitting  parameters  related  to  the  ge¬ 
ometry  of  the  porous  material.  The  form  of  this  relation  was  initially 
proposed  by  Gebart  [41]  for  the  permeability  of  ordered  two- 
dimensional  arrangement  of  cylinders. 

The  reported  values  of  (j>c,  C\,  and  C2  for  three-dimensional 
fibrous  porous  materials  with  random  structure  are  0.074,  0.491, 
and  2.310,  respectively  [27],  In  this  work,  the  fibers  of  reconstructed 
GDL  samples  have  a  preferred  direction,  thus,  permeability  is 
different  in  the  in-  and  cross-plane  directions,  where  both  are  also 
different  from  permeability  of  random  fibrous  structures. 

In  fitting  the  Eq.  (7)  to  our  numerically  predicted  data,  we  keep 
the  values  of  0C  and  C2  constant  and  equal  to  those  previously  re¬ 
ported  for  random  fibrous  structure.  New  values  of  Ci  for  in-  and 
cross-plane  directions  for  three  different  values  of  fiber  volume 
percentage  are  predicted  from  curve  fitting,  and  are  presented  in 
Table  2.  As  the  behavior  of  samples  at  the  high  porosity  end  is 
erratic,  we  only  include  samples  of  porosities  smaller  than  0.85  in 
the  curve  fitting  process. 

Fig.  7  shows  the  effect  of  fiber  percentage  on  the  predicted 
values  of  in-  and  cross-plane  tortuosity.  Decreasing  the  fiber  per¬ 
centage  (adding  more  binder  &  PTFE)  increases  the  predicted  values 
of  tortuosity.  The  effect  of  fiber  percentage  on  tortuosity  values  is 
more  intensified  at  low  porosities,  and  is  almost  negligible  for 
porosities  larger  than  0.80. 

The  increase  in  the  tortuosity  by  increasing  the  volume  per¬ 
centage  of  binder  &  PTFE  (decreasing  the  fiber  volume  percentage) 
can  also  be  explained  by  considering  the  shape  of  fibers  and  binder 
&  PTFE  material.  Large  bulk  volumes  of  binder  &  PTFE  material 
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Fig.  7.  Predicted  values  of  in-plane  (top)  and  cross-plane  (bottom)  tortuosity  of 
reconstructed  GDL-like  samples  with  heterogeneous  porosity  distribution  and  three 
different  fiber  percentages  of  50%,  75%,  and  100%. 
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effect  in  samples  with  two  different  fiber  volume  percentages  of  50 
and  75%.  Fig.  8  shows  the  predicted  values  of  in-  and  cross-plane 
permeability  for  samples  with  homogeneous  and  heterogeneous 
distributions  of  porosity  and  two  different  fiber  volume  percent¬ 
ages.  It  can  be  seen  that  for  both  in-  and  cross-plane  directions,  and 
for  both  volume  percentages  of  fibers,  the  permeability  of  hetero¬ 
geneous  samples  is  consistently  lower  than  homogeneous  samples 
for  the  entire  range  of  the  examined  porosity. 

4.5.  Comparison  with  experimental  measurements 

In  order  to  validate  the  methodology  and  have  a  comparison 
with  previously  reported  values  of  permeability  of  different  GDL 
samples,  we  compared  the  numerically  predicted  in-  and  cross¬ 
plane  permeability  values  of  samples  with  three  different  fiber 
volume  fractions  against  some  of  the  available  experimental  data  in 
the  literature  (Fig.  9)  [14-16,19,38,42-45],  As  some  of  the  experi¬ 
mental  works  have  not  reported  the  radius  of  fibers  in  their 
experimental  samples,  we  have  reported  the  absolute  values  of 
permeability.  It  can  be  seen  that  for  both  in-  and  cross-plane  di¬ 
rections,  the  numerically  predicted  values  are  in  very  good  agree¬ 
ment  with  the  experimentally  measured  values.  One  should  note 
the  scatter  in  the  experimentally  reported  data  and  the  linear  scale 
of  the  vertical  axis. 

5.  Conclusions 
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Fig.  8.  Effect  of  porosity  heterogeneity  on  the  predicted  values  of  normalized  in-plane 
(top)  and  cross-plane  (bottom)  permeability  of  reconstructed  GDL  samples  with  het¬ 
erogeneous  and  homogeneous  porosity  distributions,  and  two  different  fiber  per¬ 
centages  of  50%  and  75%. 


distort  the  flow  pattern  and  streamlines  more  than  the  slim  fibers, 
thus,  increase  the  tortuosity  of  fluid  motion. 

4.4.  Effect  of  porosity  heterogeneity  in  samples  with  binder  &  PTFE 

In  Section  4.2,  we  studied  the  effect  of  porosity  heterogeneity  on 
samples  consisting  of  only  fibers.  In  this  section,  we  study  this 


We  used  a  pore-level  numerical  modeling  approach  to  investi¬ 
gate  the  effect  of  porosity,  porosity  distribution,  and  the  volume 
percentage  of  binder  8r  PTFE  material  on  the  bulk  hydrodynamic 
properties  of  reconstructed  GDL-like  samples.  It  was  shown  that 
the  effect  of  porosity  heterogeneity  is  negligible  in  comparison  to 
the  other  two  factors  for  all  the  samples  studied  here,  which  are 
porosity  and  fiber  volume  percentage. 

Increasing  the  volume  percentage  of  binder  &  PTFE  in  the 
structure  (reducing  the  fiber  percentage)  increases  both  in-  and 
cross-plane  permeability  and  tortuosity  values.  A  previously  re¬ 
ported  relation  for  the  permeability  of  fibrous  porous  materials  was 
adopted  and  fitted  to  the  predicted  values  of  in-  and  cross-plane 
permeability  for  three  different  fiber  volume  percentages.  The 
magnitude  of  a  fitting  parameter  in  that  relation  was  reported  as  a 
function  of  fiber  volume  percentage.  The  predicted  values  of  in-  and 
cross-plane  permeability  show  a  very  good  agreement  with  the 
available  experimental  data  in  the  literature. 
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Fig.  9.  Comparison  of  the  numerically  predicted  values  of  (top)  in-plane  (bottom)  cross-plane  permeability  of  sample  with  heterogeneous  porosity  distribi 
fiber  percentages  of  50%,  75%,  and  100%  against  some  of  the  available  experimental  data  in  the  literature. 


three  different 
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